Two-dimensional frustrated Heisenberg model: Variational study. 
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The stability of the ferromagnetic phase of the 2D quantum spin-| model with nearest-neighbor 
ferro- and next-nearest neighbor antiferromagnetic interactions is studied. It turns out that values 
^\ ' of exchange integrals at which the ferromagnetic state becomes unstable with respect to a creation 

of one and two magnon are different. This difference shows that the classical approximation is 
ON 1 inapplicable to the study of the transition from the ferromagnetic to the singlet state in contrast 

with ID case. This problem is investigated using a variational function of new type. It is based on the 
boson representation of spin operators which is different from the Holstein-Primakoff approximation. 
This allows us to obtain the accurate estimate of the transition point and to study the character of 
| the phase transition. 
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In recent years the investigation of two-dimensional frustrated Heisenberg model is of great interest. This is mainly 
caused by studies of magnetic properties of superconducting cuprates. The so called J\ — J 2 — J3 model was studied 
by different methods in the case of completely antiferromagnetic interactions J\, J 2 , J3 > [1-8]. An existence of 
disordered phases at the T — and nontrivial ground states in this papers is generally assumed. 

The model with ferromagnetic interactions of nearest neighbors and antiferromagnetic interactions of next nearest 
neighbor spins has been investigated much less. The Hamiltonian of this model has a form: 

ff = -]T(S n .S n+a -i) + J 5> n -S n+d -i) ; (1) 

c • 
o . 

where vectors a and d connects nearest sites and nearest sites along the diagonal line respectively, n - a site number 
and J > 0. The model (1) is a special case of the J\ — J 2 — J3 model with J\ = — 1, J 2 = J and J3 = (J3 - is the 
exchange integral of the next nearest neighbors along X and Y axes) . This model is the simplest example of the 2D 
frustrated system. 

The ground state of the Hamiltonian (1) is ferromagnetic at small J and the energy E = 0. But this state becomes 
ON , unstable at some critical point J c . In the classical approximation J c is equal to 1/2 and the ground state at J > J c 
corresponds to two independent sublattices with Neel order. However, quantum fluctuations can change this situation. 
This question will be discussed in this paper. 

To clarify the situation let us compare the 2D model (1) with its ID version which has been studied in detail 
' previously [9-11]. In the ID model the transition from the ferromagnetic ground state to the spiral one takes place 
at J c = 1/4. In a recent work by two of present authors [11] the character of this transition has been investigated 
by the perturbation theory in the small parameter (J — J c ) and the classical state has been used as zeroth-order 
approximation. The ground state in all cases turned out to be either with S = or with S = S max , and the transition 
1 ■ from the ferromagnetic state to the state with S = occurs by passing the states with intermediate spins. Quantum 
fluctuations do not change the critical point J c , which coincides with its classical value. 

This method has been applied also to the study of the 2D model with Jj = — 1 and J 2 , J3 > in the vicinity of 
the phase boundary 2J 2 + 4J3 = 1 (Fig.l), which determines the region of stability of the ferromagnetic state in the 
classical approximation. At J 2 < 0.36 the situation is similar to the ID case [12]. However, at J 2 — > 0.36 the second 
order of the perturbation theory in the small parameter describing a deviation from the phase boundary diverges. 
It proves that quantum fluctuations change the classical phase diagram itself. The energies of one- and two-magnon 
states are shown in Fig. 2 as functions of the parameter J for model (1) (J 2 = J and J3 = ). The critical values 
Jc{Smax — 1) and JciSmax — 2) of the instability of the ferromagnetic state with respect to a creation of one and two 
magnons are equal to J c (S max — 1) = 1/2 and J c (S max — 2) = 0.4082 respectively (two-magnon state energy was 
found by the numerical solution of the corresponding Schrodinger equation). The difference between J c (S max — 1) 
and J c (S max — 2) shows that the classical approach is inapplicable to the study of the stability of the ferromagnetic 
state. In this respect the situation is essentially different from those in the ID or in the 2D cases at J 2 < 0.36, when 
the critical point J c is independent on the number of magnons and coincides with its classical value. It is natural 
to assume, that the critical value J c {S m nx — 3) for model (1) is less than J^iSmax — 2), and the critical value J c (0) 
corresponding to the state with S — is a true point of the phase transition. So, the classical approximation has 
proved to be unsuitable for the study of the phase transition in the model (1) and for the determination of the critical 
value J c and, therefore, another approach is needed. 
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FIG. 1. T = phase diagram of the Ji, Ji, Jz model with Ji = —1. The bold line is a boundary between the ferromagnetic 
and singlet phases in the classical approximation, the bold+thin lines - in present approach. 
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FIG. 2. Energies of one- and two-magnon states. 



One can obtain a crude estimate of the energy and the critical point J c (0) by using a product of the ground state 
wave functions of two independent antiferromagnctic sublattices as a variational wave function (VWF). This VWF 
gives for the critical point J c (0) a value: 

J c = y^Ye - °- 428 ' ( 2 ) 

where e is the ground state energy per site of the 2D antiferromagnetic Heisenberg model, and we use the most 
accurate numerical estimations of e = —0.669 [13,14]. Of course, this approach is too poor, and the obtained value 
J c (0) is greater than J c {S max _ — 2). 

In the present work we propose a new type of VWF, which allows us to get more accurate estimations of J c and 
to study the character of the phase transition. This approach is based on a boson representation of the Hamiltonian 
(1) which is different from the Holstcin-Primakoff one. In contrast with the spin wave theory (SWT) the proposed 
method is variational. 

This article is organized as follows. In the next section we demonstrate the method by the application it to the 2D 
Heisenberg antiferromagnetic model. In section 3 the stability of the ferromagnetic phase of the frustrated model is 
studied and the data are discussed. 



2 



II. THE METHOD 



To illustrate the main features of our approach, let us consider the 2D HAF model: 

# = ^S n -S n+a , (3) 



n,a 



where n is a site number of the 2D lattice and a vector a connects nearest sites. 

It is convenient to rotate the local coordinate system of one of sublattices by the angle ir in XZ plane: 



H — ( ^" ' ^n+a + ' ^n+a ' ^n+a) (4) 



The transformation from spin-operators to bose-ones is defined by: 



, A , A 

where bj~ are bose-operators, Ni= 6?~6j, and the operator function 9(N) is: 

It is evident that this transformation preserves all commutation relations for the spin operators. The states with 
different numbers of bosons on each site are effectively separated into equivalent unconnected pairs: 



{ 



S z \2m) = -\ , S+ 1 2m) = |2m+ 1) , S~ |2m) = 
5 z |2m+l) = i , S+|2m+l)=0 , S 1 " |2m + 1) = |2m) 



(7) 



As a result of the transformation (5) the Hamiltonian (4) takes the form: 

V V 7V„ V iVn+a 



n,a 



This Hamiltonian, as well as the original one, can not be solved exactly. As a trial wave function for Hamiltonian 
(8) we choose : 



|*>=exp^X>(i-j) 6+6+ j 



Ob) , (9) 



where the function A(i — j) will be found by the minimization of the total energy. 

The vacuum state in Eq.(9) corresponds to the state of the Hamiltonian (4) with all spins pointing down (or 
corresponds to a "chess" arrangement of spins for the Hamiltonian (3)). 

To calculate the ground state energy we need to calculate expectation values of all terms in the Hamiltonian (8) with 
respect to VWF (9). At first we calculate the terms corresponding to the spin interactions along the horizontal line 
(X axe). Owing to the translational symmetry of VWF (9), all of these terms give equal contributions to the energy, 
and, therefore, it is sufficient to calculate terms corresponding to the interaction Si • S2 in the original Hamiltonian 
(3), where 1 and 2 are nearest neighbors along the horizontal line. 

At first we represent all factors } in the Hamiltonian (8) in the form: 

V Nj 

OO 

' =4= / e-^^dtj (10) 



3 



Then the expectation value of the second term of the Hamiltonian (8), corresponding to the term —SfS^ + SfS^ 
in the Hamiltonian (4), takes the form: 



TJ X V 

H 12 



dtidt 2 
8tt 



(11) 



Thus, we need to calculate expectation values of the type: 



and 



where r\ and can be: 



A\2 = 



r - — t- 



Tj = tj + in 



.7 = 1,2 



(12) 



(13) 



(14) 



In order to decompose the boson quadratic form in Eq.(9) we use the Hubbard-Stratanovich transformation. Now, 
the VWF (9) takes the form: 



i*>=t^ / n<^ EVwj ^ io b) 



where N is a number of lattice sites. 

Now we can calculate the expectation values (12,13): 



Cl2 



-OO 1 

OO ' 

— oo 1 



(15) 



(16) 



where we use the notations: 



-oo 1 

' oo • 

/ Ild&drnWtr, 

— oo 1 



Wt„ = exp A ij 1 (6£i + ViVj) + ^ j ' 



Pj = l~e~ r \ 



J = 1,2 



Using the identity: 



Eqs.(16,17) can be written as: 



66 +vm = 



dp 3 



a P3(?l?2+r)1772) 



P3=0 



C12 = ( e 



-2(^ir)i+^2J)2) 



(17) 



(18) 



(19) 



(20) 
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where matrix has the form: 



Bij - 



( p' 2 + a+ -p' 3 -vt-2 \ 

p' 2 + a- -p' 3 + 

\ -P3 + V1-2 p[+a- ) 



Substituting Eqs.(30,31) into (21), we find: 

d\2 = e~ ri ~ r2 (l4- 2 fl(Pl,P2) f-(Pl,P2) - Ml~-2 f+(Pl,P2) f-(Pl,P2)) , 

where 

f±{PUP2) = [{l+p 1 <7 ± ){l+ P 2<? ± )-piP2(}lt-2?Y k 

There are four terms in (11), corresponding to the different values n, r 2 and p\, p 2 in (14) and (19). 



p a = 1 + ae *i 
9/5 = 1 + /3e-*= 



a,/3 = ±1 



(31) 



(32) 



(33) 



(34) 



Using the notation (34), Eq.(32) and the definition of functions / + and /_ in (33), Eq.(ll) can be written as: 



&12 - + h 12i 



(35) 



(36) 



a,/3=-l 



oc 

_ ff dhdt 



2 -t 2 -t 2 

e ' 



- ^ f+(p a ,qp) ft(p a ,qp) 

a,0=-l 



(37) 



We note, that the expectation value C12 can be obtained from Eq.(30) by setting pi = p 2 = 2 and p 3 = 0. Hence, 
we obtain: 



£?2 = Jci 2 = i/ + (2,2)/_(2,2) 



(38) 



The contribution to the energy of the terms corresponding to nearest neighbors along the horizontal line in the 
Hamiltonian (8) is: 



E12 — —Ef 2 + E 12 — E\ 



(39) 



One can verify by rotating the local coordinate systems and repeating the calculations (5-30) for the transformed 
Hamiltonian, that the values Ef 2 and E\ 2 correspond to expectation values : 



E?2 = (S?S5) 



E\ 2 = (SfS y 2 ) 



(40) 



The calculation of the expectation values of the terms of the Hamiltonian (8) corresponding to the nearest neighbor 
interactions along the vertical line (Y axe) Si • S3 can be carried out in the analogy with (9-30). In this case, we 
must only substitute the functions fJ,f_ 3 for fif_ 2 into Eqs. (36-38). Then we obtain the contribution to energy of the 
vertical neighbor terms: 



£13 - --E13 + Ef 3 - E* 3 
The total energy of Hamiltonian (8) can be written as: 

E = E\2 + ^13 



(41) 
(42) 
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Thus, we find energy as a function of parameters ct ± , /if_ 2 , /if_ 3 . These parameters are functions of w(k) [or of 
A(i — j) in VWF (9)]. Now we need to minimize the energy with respect to w(k): 

dE ai + bi cos k x + c\ cos k y 0,2 + &2 cos k x + C2 cos k y 

a^ = K - 1) 2 + (^k + 1) 2 = ' ( } 

where 

n — dE h — dE n — dE 

ai "^' 5l ~^rv Cl "^? 

n — JL§L h — 9E o — 9E 

We can define the functional form of w(k) from Eq.(43): 
where 



. / 1 + a 2 cos fcj; + a 3 cos k v 

g(k) = ai J — — f-, (45) 

y l + «4 cos k x + CX5 cos k y 

and ai (i = 1, 5) are variational parameters. 

Thus, we reduce the problem of the energy minimization over the variational function oj(k) to the energy minimiza- 
tion with respect to five variational parameters ai . This procedure has been performed numerically. It gives the final 
result for 2D HAF: 

ai = 1, a 2 = a 3 = -a 4 = -a 5 = -0.445, E = -0.640 (46) 
Besides, both terms in Eq.(42) give equal contributions to the energy: 

E12 = E13 

Obtained ground state energy for the 2D HAF (3) is 4-4,5% higher than the most accurate results obtained by 
various numerical techniques [13,14] and is in a good agreement with energies obtained by different VWF methods 
[15-17]. 



III. RESULTS AND DISCUSSION 



Now we apply the proposed approach to the frustrated model (1). The ground state of the Hamiltonian (1) is 
ferromagnetic at small J. In the classical approximation two-sublattice Neel state is realized at J > |. The energy 
of this state is: 

if=-0 < 47 > 

A spiral state (the incommensurate phase with the momentum (Q,Q), cosQ = j-j) has higher energy at J > |: 

-U-\f (48) 

although it tends to zero at J — > ^ as well. 

Energies of other phases are considerably higher. The energy of the dimer phase, for example, equals to = | at 
J = i. Hence, it is natural to consider Neel two-sublattice and the spiral states as main candidates for the ground 
state at J > J c . We will calculate quantum corrections to the classical energy for these states using the VWF (9). 
Before that, we make the following remark. 

The VWF (9) has no rotational symmetry ( it is not the eigenfunction of S 2 ), and, therefore, the ground state 
energy calculated in this approximation depends on the choice of the local coordinate system for the Hamiltonian 
(1). In other words, if we rotate spin operators S — > q and perform a transformation to bose-operators according to 
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FIG. 3. The dependence J c on the angle tp. 



Eq.(5), then, generally speaking, obtained energies will be different in spite of the equivalence of the Hamiltonians H s 
and Hence, parameters of the rotation can be considered as variational ones. 

For the consideration of the two-sublattice Ncel phase let us rotate the coordinate system so that there are the 
angle ir in XZ plane between two nearest neighbors on each sublattice (as in the 2D HAF case) and the angle ip 
between sublattices in XZ plane. We will keep the angle ip as a variational parameter (we note, that the energy is 
infinitely degenerated with respect to (pin. the classical approximation). In this case, the Hamiltonian (1) takes the 
form: 

n 

+ Sinp (££ +a- - C„^ +ax - «S«S+a, + <^n+aj + &Z+a. + *n+aj 

+ J ^2 (~^n+d + ^n^n+d - ^n^n+d _ J J i 
n,d ^ 



(49) 



where vector d connects nearest sites on diagonal line. 

Expectation values of the terms in the Hamiltonian (49), such as (<r£ • <^), • and (<r„ • Cm) are determined 
by Eqs. (36,37,38) respectively, where /i^_ m and /u~_ m are defined by Eq.(29). 

The VWF (9) contains only terms that involve even numbers of boson operators and, therefore, all expectation 
values such as (<;* • ? m ) in (49) are equal to zero: 



In this case there are nine variational parameters in Eq.(45): 



g(k) = a M 



1 + a.2 cos k x + a3 cos k z + ot4 cos(k x + fcz) + «5 cos(fca; — k z ) 

1 + OiQ COS fcj; + Qf7 COS fc z + «8 COs(fer + k z ) + Ctg COs(k x — k z ) 



(50) 



(51) 



The critical value J c is defined by the condition that the ground state energy is negative for J > J c . To find J c we 
minimize the following ratio: 



Jc = 



(52) 



where E\ and E 2 are the contributions to the energy from the first and the second terms in the Hamiltonian (49): 

E=-E 1 + JE 2 (53) 
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FIG. 4. Energies of the collinear (solid line) and spiral (dashed) states. Dotted line is the SWT energy of the collinear phase. 

The dependence of J c on <p obtained by a minimization with respect to parameters ai is shown on Fig. 3. As it can 
be seen from Fig. 3, the minima of J c are reached at ip — or (p — it. The corresponding states belong to the so called 
collinear phase. At ip = we have: 

ai = 1, a 2 = a 6 = —0.7, 0:3 = —017 = 0.64, 

a 4 = a 5 = -a 8 = -a 9 = -0.36, J c = 0.4056 (54) 

The calculation of the energy with use of the VWF (9) for the spiral phase can be carried out in a similar way. 
It turns out that quantum fluctuations do not shift the transition point and change the coefficient in the quadratic 
dependence in (48) only. The dependences of the energies of the collinear and spiral states are shown on Fig. 4. For 
comparison, we show the spin wave theory results as well. The SWT energy is not a variational one and is defined 
up to J = 0.52, where the sublattice magnetization vanishes in the SWT approximation. 

Thus, the singlet collinear phase is the ground state of the model (1) for J > J c . 

A characteristic feature of the considered model (1) is the fact, that the critical value J c , (the point of the instability 
of the ferromagnetic state), depends on S and J C (S) < J C (S + I). In this respect this model differs from the ID case 
where J c does not depend on S. Such ID behaviour takes place also in a more general 2D J2 — J3 model in the vicinity 
of the classical phase boundary and at J2 < 0.36. In this case, quantum effects do not change the classical boundary 
of the instability of the ferromagnetic state. The region near the boundary can be considered in the framework of 
the perturbation theory, and the transition from the ferromagnetic state to the spiral one occurs in this case. The 
character of the transition changes at J2 > 0.36, and the critical parameter J c , corresponding to E(S) — 0, depends 
on S. In this case the transition from the ferromagnetic to the collinear state is realized. The boundary of the stability 
of the ferromagnetic phase has a form shown on Fig.l. 

In conclusion, we have studied the transition from the ferromagnetic to the singlet state in the 2D frustrated spin 
model. The transition region is characterized by strong quantum fluctuations and can not be described by the classical 
approximation. To study the behaviour of the system close to the ferromagnetic boundary we have proposed new 
approach based on the bozonization of the spin operators. This approach is different from the Holstein-Primakoff 
method and is variational. We believe that the proposed method can be used to the study of other Heisenberg models 
with the frustration. 
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